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Abstract 

In this work, we analyze the Born, Bogoliubov, Green, Kirkwood and Yvon (BBGKY) hierarchy 
of equations for describing the full time-evolution of a many-body fermionic system in terms of its 
reduced density matrices (at all orders). We provide an exhaustive study of the challenges and 
open problems linked to the truncation of such hierarchy of equations to make them practically 
applicable. We restrict our analysis to the coupled evolution of the one- and two-body reduced 
density matrices, where higher order correlation effects are embodied into the approximation used 
to close the equations. We prove that within this approach, the number of electrons and total energy 
are conserved, regardless of the employed approximation. Further, we demonstrate that although 
most of the truncation schemes available in the literature give acceptable ground state energy, when 
applied to describe driven electron dynamics exhibit undesirable and unphysical behavior, e.g., 
violation and even divergence in local electronic density, both in weakly- and strongly-correlated 
regimes. We illustrate and analyze these problems within the few-site Hubbard model. The model 
can be solved exactly and provides a unique reference for our detailed study of electron dynamics for 
different values of interaction, different initial conditions, and large set of approximations considered 
here. Moreover, we study the role of compatibility between two hierarchical equations, and positive- 
semidefiniteness of reduced density matrices in instability of electron dynamics. We show that even 
if the used approximation holds the compatibility, electron dynamics can still diverge when positive- 
definitiveness is violated. We propose some partial solutions of such problem and point the main 
paths for future work in order to make this approach applicable for the description of the correlated 
electron dynamics in complex systems. 



* These authors have contributed equally in this work. 
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I. INTRODUCTION 



ectron 

3, Is). 



The description of dynamical processes in many-electron systems brought out of equi- 
librium requires a proper description of static and dynamical correlation effects. Theoret- 
ical study of many fundamental processes of interest such as attosecond dynamics, high- 
intensity laser phenomena, light-induced phase transitions, harmonic generations, etc., re- 
quire a proper theoretical framework which are both accurate and tractable. 

Time-dependent Hartree-Fock (TDHF) theory [l-3] was among the first methods devised 
to simulate dynamics of many-body systems. More advanced methods such as Kadanoff- 
Baym equations 0, 0] and the Keldysh technique ||3] for the Green's function enable us to 
study non-equilibrium phenomena more accurately, once an approximation for the e 
self-energy is chosen properly. They are, however, computationally very demanding 

Alternatively, time-dependent density functional theory (TDDFT) 0, [lo| is becoming a 
popular option for its performance, accuracy and potential scalability. It describes electron 
dynamics in terms of the electronic density by mapping an interacting system to an auxiliary 
noninteracting system (Kohn-Sham system) with an effective potential that produces the 
same time-dependent density. This is ideal for massively parallel implementations as in 
its time-dependent Kohn-Sham (KS) formulation the evolution of each of the KS orbitals 
is nearly independent of the others. However, lack of proper exchange-correlation (xc) 
functionals for time-dependent systems hamper its applicability. Most of the time-dependent 
approximations that are used so far, are adiabatic extensions of DFT approximations that 
disregards the non-locality and memory dependence of the time-dependent xc potential. 

In this situation, reduced density matrix (RDM) theories offer a promising framework 
to deal with time-dependent phenomena. The equation of motion for each RDM can be 
straightforwardly derived from the time-dependent Schrodinger equation which contains the 
corresponding and one order higher RDM. The whole set of these interrelated equations form 
the so-called BBGKY hierarchy since a basically similar hierarchy was initially invented and 



llMl4j in classical statistical 



developed by Born, Bogoliubov, Green, Kirkwood and Yvon 
mechanics. 

As it is not possible to confront the whole hierarchy, one must truncate it at some level. 
For instance, if one approximates two-body RDM in terms of one-body RDM in the first 
equation, one arrives at the time-dependent version of reduced density matrix functional 



theory (TD-RDMFT). Similar to the TDDFT, most of the approximations used in TD- 
RDMFT are adiabatic extensions of the existing ground state ones 15 and even though 



IS, 



19|, 



they can successfully describe the ground state of some strongly correlated systems 
they suffer from flaws such as lack of memory when we apply them in time domain. As an 
example, assuming that these approximations have even the right form, the sign adopted in 
each term generally becomes a time- dependent phase when we extend them over the time 
domain. Ignoring this fact by using fixed signs may cause problems such as time-independent 
occupation numbers 20|. Furthermore, majority of these approximations do not necessarily 
conserve total energy of a system. 



Some of these deficiencies will be cured if we consider propagating the first two equa- 
tions of the hierarchy by approximating the three-body RDM. However, this prove to be 



a nontrivial task and in fact there are earlier attempts in nuclear dynamics [21|, |22] which 
show that these coupled equations can violate the inequalities related to probabilistic in- 
terpretation of RDMs for fermions. For example, the eigenvalues of one-body RDM must 
remain, in principle, between zero and one 23] but in practice they observed violation of 
these bounds, showing the non-fermionic nature of the corresponding RDM. Such behaviors 
were unexpected and it was claimed to be related to the violation of the relations between 
different orders of reduced density matrices. 



In this paper, we study in detail the performance of such approach for different truncation 
schemes and show that the truncated set of equations may lead to instability and in many 
cases even divergence (in electronic density, occupation numbers, etc.). We mention the 
specific properties of approximations that are responsible for these unphysical results. We 
will show that lack of properties such as positive-semidefiniteness also plays a crucial role 
in this failure. In addition, this study prompts one to be aware of the same issues which 
may arise in building approximations in TD-RDMFT. The article is divided into three 
sections. In next section, we give the theoretical background for the BBGKY hierarchy and 
different approximations to truncate it. In Section II III we present the results and analyze 
the phenomenon using different approximations, and finally we conclude the work in Section 
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II. THEORETICAL BACKGROUND 



A. N-Particle System and Reduced Density Matrix: Definitions and properties 

First in order to have a compact notation, we introduce two collections of space-spin 
coordinates as 

X n = (xi . . . x n ) ; X n = (x n+ i . . . xjy). (1) 

In this notation, Q(XN,t) denotes the normalized wave function of the system. Here, we 
note that throughout this work we employ atomic units. 

Now, we consider a system with N identical particles described by the time-dependent 
Hamiltonian 

N 1 N 

= 2 E^- ( 2 ) 

8=1 i^j 

The last term describes the two-particle interactions Uij = ?7(xjXj) where x includes both 
space coordinates, r, and spin coordinates, a, of the particles. Usually U is spin-independent 
and has the form [/(xx') = w(\r — r'|). The one-body part, hi = /i(x,j,t), will be time- 
dependent and of the form 

h(x,t) = -^V 2 + ^(x,t) (3) 

where v is a general time-dependent external field. We can define the n-body reduced density 
matrix, Y^ n \ of such system as 

AT' /* 

T^\X n ,X' n ,t) = — J dX n $(X n ,X n ,t)&(X' n ,X n ,t) (4) 

where dX n = dx n+ i . . . cZxjv and J dx. — J2 a I 

Based on the above definition, several important properties of RDMs follow. First, dif- 
ferent levels of RDMs are connected to each other by 

J ^r^fl^.Xvi.t) = (N-n)T^ n \X n ,X' n ,t) (5) 

which we refer to it as partial trace relation. Consequently, if r^ n ^ is available, all RDMs 
with lower order can be calculated straightforwardly. 

Another important feature is that Eq.flJ]) implies all RDMs are positive-semidefinite which 
refers to the fact that all eigenvalues of RDMs are always equal to or greater than zero. For a 
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given order RDM, these eigenvalues and their corresponding eigenvectors can be calculated 

as 

J dX' n r^(X n ,X' n ,t) gi (X' n ,t) = Ai(t)(fc(X nj *), Xi{t) > 0. (6) 

Conventionally, the eigenvectors and eigenvalues of are called natural orbitals and nat- 
ural orbital occupation numbers, and of are called geminals and geminal occupation 
numbers, respectively. 



Moreover in the case of fermionic particles, the Pauli exc 
orbital occupation numbers to be less than or equal to one 



usion principle enforces natural 



231 ] . Thus, they have to remain 



between zero and one. This corresponds to the fermionic inequality. 



B. Equation of Motion, BBGKY Hierarchy and Conservation Laws 

Provided the initial state of the system is given, its time evolution is completely described 
by the time- dependent Schrodinger equation (TDSE). Nonetheless, the solution of the TDSE 
is not usually possible, except for systems with few particles, due to the many degrees of 
freedom of the system. In contrast, most quantities of interest are n-body observables which 
can be obtained from the n-body reduced density matrix, T^ n \ and in most systems of 
relevance it holds that n <C N. It is therefore natural to derive the equations of motion for 
reduced density matrices. Thus, using Eq.fllJ) together with the TDSE, we get 

(i d t - + H v ...n>) r (n) (x n , x' n , t) = 

n „ 

/ dx n+1 (f/( Xi x n+1 ) - f/(x:x n+1 ))r(" +1 )(x nXn+1 ,xx+i^) (7) 

i=l J 
where we defined 

n 1 n 

= + (8) 

The entire set of iV equations for reduced density matrices forms the well-known BBGKY 
hierarchy. The explicit form of the first two equations are 

(i d t -h 1 + h v ) x' l7 1) = J dx 2 (U (xix 2 ) - U (x / lX 2))r( X ix 2 , y^x 2 , t) (9) 

and 

(i d t - H 12 + H V2 <) r(x 1 x 2 , x' lX2 , t) = 

J dx 3 ([/( Xl X 3 ) + f/(x 2 X 3 ) - [/(x^Xs) - f/( X2X3 ))r (3) (x lX2X 3, x'xX^Xs, *) (10) 
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where 7 = I^ 1 ) and T = T^. As it is customary in the literature 24[, we call the right- 
hand side of Eq.f llOp the three-body collision integral and use S to refer to it. In general, 
each equation of the hierarchy relates a given- order of RDM, T^, to one order higher 
RDM, r( n+1 \ In order to make the BBGKY hierarchy practical, we must truncate it at 
some level, n, by reconstructing the r^ n+1 ^ as a functional of lower-order RDMs. Although 
such reconstruction in terms of the one-body time-dependent density matrix is, in principle, 
conceivable by virtue of the Runge-Gross theorem 0, [25}], there is no practical method 
available to find the exact functional and we have to use different approximations. 

In this work, however, we will only propagate 7 and T since they are sufficient to calculate 
the dynamics of all one- and two-body observables. For instance for the case of total energy, 
the expectation value of the one- and two-body part of the Hamiltonian, Ei(t) and E 2 (t), 
can be written as 

N 

Ei(t) = X>*> = / rfx'/i(x',t)7(x,x',t)| x=x , (11) 

and 

1 N 1 r 

E2(t) -B^) = 2 dX 2 U(X 2 )F(X 2 ,X 2 ,t). (12) 

i,3 

To this end, we need to truncate Eq.f llOp by approximating I^ 3 ) in terms of 7 and Y. Several 
of these approximations will be discussed shortly. 

At this point we highlight an important property of BBGKY hierarchy and the effect 
of truncation on it. As a direct outcome of Eq.(j5]), different levels of the hierarchy are 
compatible; namely, equations in the higher levels of the hierarchy are reducible to the 
lower-level ones. We refer to this link between equations as compatibility condition that 
preferably should be fulfilled by a good approximation. Thus, compatibility signifies that 
the highest equation is equivalent to the whole BBGKY hierarchy. This is not surprising 
since the highest equation is basically the original Schrodinger equation. However, when we 
truncate the hierarchy by introducing an approximation for the partial trace relation 
between I^ 3 ) and V does not necessarily hold and thus it generally breaks the compatibility 
between Eqs.Q and ffTU]) (see subsection IIII CI for some exceptions). Consequently, when 
we truncate the BBGKY hierarchy, we have two generally distinct options to propagate the 
equations which should be equivalent if the truncation approximation satisfies compatibility: 

i) Propagating two coupled equations: We can evolve both 7 and Y by solving Eqs.([9]) 
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and (flOl) together as coupled equations since the two equations most likely are not 
compatible anymore after approximating I^ 3 ) in Eq. OlOp . 

ii) Propagating only second equation: To avoid the problem of compatibility between two 
equations, we can evolve only Eq. fllOj) . Then we assign 7 to be the partial trace of T 
and denote it as 7 r to distinguish it from general 7. It mathematically reads 

7 r (xi,x' l7 t) = - J rfX2r(x!X2,xix2,t). (13) 

In this way, we prevent the complication of dealing with two coupled equations. 

The difference between these two approaches lies in distinction of 7 and 7 r . To see this, 
we derive the equation of motion for Eq. ffTBl by using Eq. ffTUl) . We have 



i^7 r (xi,xi,t) = j^-j J d.-x 2 id t r(X 2 ,X 2 ,t)\ X2=x > 2 

i-^y dx 2 (H 12 -H V2 ,)T(X 2 ,X 2 ,t)\ 



N 
4 



X 2 =X, 



^— j" J rfX2dx 3 (f/(xiX 3 ) - f/(xiX 3 ))r^(xiX2X3,X / 1 X2X3,t) (14) 



N 

where the approximated 3-body RDM, T a p P , does not necessarily integrate to (N — 2) T as 
we pointed out earlier. Now, using the explicit form of Hi 2 , we can rewrite Eq.( fl4|) as 

Z<9 t 7 r (xi,xi,t) = ( h l ~ ^l')7r( x l 5 x 'l^) 

+ J dx 2 (tf( Xl x 2 ) -f/(x / 1 x 2 ))r(x 2 ,X2,t) 

~ j rfx 2 {v\ - v' 2 2 )r(x 2 ,x 2 ,t)| X2=x , (15) 



where 

f(x 2 ,x 2 ,t) = ^ T 

The last term in Eq.f fT5|) can be written as a total divergence as follows 



r(x 2 ,x 2 ,t) + J dx 3 rg),(x 1 x2x 3 ,x / 1 x 2 x 3 ,t) 



(16) 



(v 2 - v' 2 2 )r(x 2 ,x 2 ,t)| X2=4 = v 2 • [(v 2 - v 2 )r(x 2 ,x 2 ,t)| X2=x J (17) 

and hence vanishes after integration over x 2 . Finally Eq. (TT3]) becomes 

(id t -h + hy) 7 r (xi,x' 1 ,t) = J dx 2 (C/(xiX 2 ) - J7(x' 1 x 2 ))f (xix 2 , x' x x 2 , t). (18) 

Obviously, the equations of motion for 7 r and 7, Eq. flTS]) and Eq.(j9]) respectively, will be 
equivalent if T — T. This derivation shows that Eqs.flH]) and ( ITS"]) , and therefore the two 
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above-mentioned approaches, are generally different as a consequence of the fact that Ta PP 
does not necessarily integrate to (N — 2) T. 

It is a priori not clear which of the two approaches is preferable. For that, let us examine 
if the approximate equations maintain important properties such as particle number and 
energy conservation. 



1. Particle number conservation 



First, we discuss the particle number conservation. To do this, we note that the diagonal 
of one-body RDM, 7(x, x, t), gives the particle density, n(x, t). Now, particle number con- 
servation means that the total number of particles, N(t) = J rfxu(x,t), is independent of 
time, i.e. dtN(t) = 0. This is guaranteed once the continuity equation holds. In fact, this is 
the case for both approaches since from either Eq.([H]) or ffTSl) . we arrive to 



t n(x,t) = <9i7(x,x',t)| 



l(V 2 -V' 2 b(x,x',t)| x , = 



-^(V 2 + V • V - V • V - V' 2 ) 7 (x, x', t) | x , 

_ V . I(V-V')7(x,x',t)| x , = 
-V-j(x,t) 



(19) 



where we defined the current density, j(x, t), as the quantity inside the braces. This follows 
immediately from of the fact that the right-hand side of both Eqs.Q and (Tl8l) vanishes for 
X! = x'j independent of T or T. Thus, the total number of particles is preserved even if we 
cut the hierarchy at the first level. 



2. Energy conservation 

Now we turn to the energy conservation. The total energy of the system is the sum of 
Eqs.( fTT|) and ( fT2l) . For the second approach where we use only Eq.l fTOl) . the time derivative 
of the total energy after some algebra (see Appendix [AJ arrives at 

= J dx3tv(x,t)n(x,t) (20) 

+ 27 / rfxi(lx2Vlf7 ( XlX2 ) ' - Vi')(r(xix 2 ,x' 1 x 2 ,t) - f (xix 2 ,x / 1 x 2 ,t))|i = i'] . 



In the absence of time-dependent potentials, v, this equation gives dE/dt = 0, provided 
r = f. Nonetheless as we discussed, this is not generally valid in the second approach and it 
means the total energy is not conserved there. In contrast, when we solve both Eqs.(l5|) and 
[TU]) together, the second term of Eq. ([2"0l vanishes automatically (Eq.f lAlOl) in Appendix IA~|) 



24. 



26j and the energy remains constant. 



The energy conservation gives us a strong motivation for preferring the first approach, in 
which we propagate both equations simultaneously, over the second and this is what we do in 
the remainder of this work. However, we look into the other case only briefly in section IIIIDI 
and show that the energy can fluctuate relatively large in time for all the approximations. 



C. Hierarchy Truncation Methods: Approximating the I^ 3 ) 

In this section we discuss in detail different truncation schemes that we have evaluated 
in the present work. One systematic way of building these approximations is called cluster 
expansion which is a method of reconstructing higher-order RDMs as anti-symmetrized 



products of lower-order ones plus a residual correlation function 26M30I] . To have a compact 
notation, first we define the wedge product as the anti-symmetrized product of p- and m- 
point functions by 

a(X p ,X' p )Ab(X p ,X' p ) = 

) Yl e (") a ( Xa i • • • X « P ' X /3i • • • X /? P ) 6 ( X «p+i • • • X "iv> X /? P+ i • • • X /3iv)- ( 21 ) 

a,j3 

Here, N = p + m, a represents all permutations of the unprimed coordinates, (3 represents 
all permutations of the primed ones, and the function e(a) returns +1 when the permutation 
a contains an even number of transpositions and —1 for an odd number of transpositions 



311 ] . For instance, the wedge product of two general one-particle matrices is 

a(xi,xi) A&(x 2 ,x' 2 ) = ^|a(xi,x / 1 )6(x 2 ,x / 2 ) - a(xi, x' 2 )6(x 2 , x'J + (22) 

a(x 2 , x' 2 )fe(xi, x[) - a(x 2 , x / 1 )6(xi, x' 2 



Now, we illustrate the cluster expansion by some examples. The first term of the ex- 
pansion of r< n > has the same form as in the noninteracting-particle picture, namely, it is an 
n-dimensional determinant of 7, with 7(x i; x^-, t) placed in row i and column j. For instance, 
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for r( 2 \ the first term reads 



2 7 A 7. 



(23) 



7(xi,x' 1 ,t) 7(xi,x' 2 ,t) 
7(x 2 ,x / 1 ,t) 7(x 2 ,x' 2 ,t) 

Now, we define a two-body correlation function, A^ 2 \ as a means of the deviation of T from 
the noninteracting form such that 



r(x 2 , x' 2 , t) = 2 7 a 7 + a( 2 ) (x 2 , x' 2 , t). 



(24) 



If we, for instance, approximate T app = 27 A 7 and replace it in the first equation of the 
BBGKY hierarchy (J9j), we recover immediately the well-known TDHF equation. 

For accordingly, we use a noninteracting particle form and add anti-symmetrized 
products of 7 with the correlation function A^ 2 ^ - that partly describe the 3-body correlation 
- plus a remainder, A®, i.e. 



r( 3 )(x 3 ,x 3 ,t) 



7(x 1 ,x' 1 ,t) 7(xx,x 2 ,i) 7(x 1 ,x / 3 ,t) 

7(x 2 ,xi,t) 7(x 2 ,x 2 ,t) 7(x 2 ,x' 3 ,t) (25) 

7(x 3 ,x' 1 ,t) 7(x 3 ,x 2 ,t) 7(x 3 ,x 3 ,t) 
3 

+ ]T (-1)^ 7 ( Xi , x ;, t) M 2) (x,, x ;, t) + a< 3 ) (x 3 , x' 3J t). 

In the second term on the right-hand side, Xj denotes the pair of variables in the set (x!X 2 x 3 ) 
complementary to Xj keeping the order of the arguments fixed; the same goes for the primed 
coordinates. For example, x 2 = (X1X3). Using the wedge product notation, we can rewrite 
Eq.(l25l) as 



= 67A7A7 + 97A A (2) + A (3) = -12 7 A7A7 + 97Ar + A (3) 



(26) 



in which we replaced the = T — 27 A 7 from Eq. fTSlI) . Similarly, we can write the 
expansion for higher-order RDMs. 

The same method has been used in the Contracted Schrodinger Equation formalism 
(the hierarchical set of equations for density matrices derived from the time-independent 
Schrodinger equation) and referred to as cumulant expansion 32l-l35l|. Nakatsuji and Ya- 
suda made the expansion more grounded by deriving it using the relation between RDMs 
and Green's functions 



331 ] . Based on these, We are now ready to discuss a number of 



approximations for r®: 
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i) Three-body collision-integral-free (3b-CIF) approximation. The simplest one rises from 
the assumption of = 0, which removes the whole right-hand-side of Eq.([9]). 

ii) Three-body-noninteracting approximation (3b-NIA). This is obtained only by consider- 
ing the noninteracting term of Eq. (1261) 

r£U M = 6 7 A7A 7 . (27) 
This gives I^ 3 ) as a functional of 7. 

iii) WC approximation. We can, of course, climb to the next level and take also the second 
term of Eq. (1261) into account which leads us to 

rg> c = -12 7 A7A7 + 9 7 Ar (28) 

where the index stands for Wang and Cassing who introduced this approximation in 
1985 j^j. This properly reduces to Eq. (T27|) when we assume T = 2 7 A 7 . 

Now, if we want to go beyond these approximations, in Eq. (l26l) must be approxi- 



33 



-|35j - for instance, by using the con- 
but we are not going to deal with 



mated. In fact, there exist some approximations 
nection between RDMs and the Green's functions 
them in this paper. 

III. DISCUSSION OF THE RESULTS 



In order to test the aforementioned approximations, we need a system for which we have 
access to its exact or nearly exact solution. The Hubbard model fits very well here since 
first of all we can solve it exactly for a few sites; and secondly, since the number of single- 
particle orbitals that build the many-body Hilbert space is limited, we can retain the full 
single-particle basis set and avoid basis set truncation errors. It is worth mentioning that 
although here we illustrate only the performance of all the approximations for the Hubbard 
model, our preliminary results for small molecules also exhibit very similar issues (work in 



progress 



36|) 



Now, we consider the case of a linear finite Hubbard chain with only nearest neighbor 
interactions. The Hubbard Hamiltonian for a finite one-dimensional system in second quan- 
tization notation is 
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H = 1 ( a l+l lCT <V + a\a a i+hv) + Y U n n n n ( 29 ) 
a,i i 

where a is a spin index, i is the site index and t and U denote hopping and on-site Coulomb 
energy, respectively. Here, t is set to unity and U gets different values to simulate different 
correlation strengths. 

To study the quality of the approximations, we must go beyond two-particle systems 
since they can be treated exactly in our formalism. In this work, however, we will avoid the 
practical complications introduced by spin in odd-number-electron systems, and perform all 
our calculations in four electrons in a four-site system but that does not affect the generality 
of our results. A more detailed discussion on numerical aspects is given in Appendix |B] and 
the code is also available upon request \^\ . 

In this part, we investigate three different approximations of T^ 3 \ namely the three-body 
collision integral free, the three-body non-interacting, and the WC approximations and 
compare them with the exact and the TDHF results. With these approximations, we have 
now a closed set of equations and as any differential equation, we need an initial state of 
the system to propagate them. To study the initial state dependence of the phenomena, we 
choose two extreme regimes of initial states to perform our calculations: Far from equilibrium 
and near to equilibrium. 

At first, we choose a far from equilibrium state as our initial state since it helps us to 
show the problem more clearly. We build such an initial state by putting four electrons in 
the two leftmost sites, i.e. 

|\&o) = a i | a i J, a 2 t a 2 4.1*-*) (30) 

where 1 and 2 refer to two neighboring sites at the beginning of the chain. In Appendix 
IB] we give the matrix form of the initial state and the Eqs.([H]) and (fTUl) that we actually 
propagate. 

Here, the equations in hand are ordinary differential equations and we solve them nu- 
merically using the Runge-Kutta method. However, to ensure the accuracy and stability 
of our results, we also used more accurate time-propagation schemes such as the fourth- 
order Adams-Bashforth-Moulton method and we found no notable difference (for a detailed 



discussion of these methods see 



38|). 



The time evolution of electronic density in the leftmost site, n(l,t), is plotted in Figfj] 
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FIG. 1. Time evolution of electronic density in the leftmost site of a 4-site Hubbard model with 
(a) TDHF (b) 3b-CIF (c) 3b-NIA and, (d) WC approximations. The exact result is also given for 
comparison. Here, to, h are set to unity and Hubbard parameters are U = 0.1 and t = 1. The four 
electrons filled the two leftmost sites initially. 



for a weak on-site Coulomb energy, U = 0.1, and for (a) TDHF, (b) 3b-CIF, (c) 3b-NIA 
and, (d) WC approximation. The plots also contain the exact result for comparison. In 
a short-time scale, we can see that all three approximations improve the quality of the 
results considerably, compared to the TDHF. However, comparing with each other, the 
approximations do not exhibit large differences. 

Figure |2] shows essentially the same results for a longer propagation time. It also shows 
how the highest and lowest geminal occupation numbers, X ma x and A m j n , behave in time. 
For the 3b-CIF in panel (a) we can see unphysical behavior around t ps 240 a.u, where the 
density acquires negative values or rises beyond two electrons in a site. The problem is more 
serious for the two other approximations since for longer propagation times, the electronic 
density starts to oscillate with amplitudes much beyond physically allowed boundaries, and 
eventually diverges as is shown in Fig. 2 (b and c). The divergence time depends on the 
correlation strength, namely on the value of U in our model, and it decreases almost ex- 
ponentially with increasing U. For example, for WC approximation, the divergence time 
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FIG. 2. Time evolution of electronic density in the leftmost site of a 4-site Hubbard model in 
a longer time scale for (a) 3b-CIF (b) 3b-NIA, (c) WC approximations. Blue lines show the 
highest and lowest geminal occupation number in time. Here, m, h are set to unity and Hubbard 
parameters are f7 = 0.1 and t = 1. The four electrons filled the two leftmost sites initially. 



changes from t ~ 532 a.u for U — 0.1 to t « 3 a.u for U = 10. It is important to note that in 
3b-NIA and WC approximations, \ max and X min start to diverge much earlier, although we 
can not immediately see the effect in neither natural orbital occupation numbers nor on-site 
electronic densities. 

Nonetheless, it is well-known that the time-evolution of a far from equilibrium state is 
generally very difficult to handle with any approximation, and particularly with the ground- 
state-tuned ones; hence, we change the initial states to be closer to the system's ground state 
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in order to investigate the generality of this phenomenon. First, we start the simulation with 
the initial 7 and T extracted from the ground state of i) the exact solution and ii) the Hartree- 
Fock approximation; and let it propagate with all three different approximations. Although 
in these cases the electronic density for the 3b-CIF does not violate physical bounds, we still 
see the divergence for other two approximations. 

Moreover, we used the method introduced by Mazziotti 39( to find the ground state asso- 



ciated with 3b-NIA and WC approximations and then used it as the initial state. However, 
since the method 39| is not totally convergent, the result is not a truly stationary state 



and even starting from such state does not bring stability to the equations and di verg ence 
appears again. We will discuss this issue in more details in our forthcoming paper 



These tests show that the divergence problem is independent of the initial state and has 
to do with the nature of the approximated equations. It is worth emphasizing again that in 
all of these approximations the continuity condition, Eq. (ll9p . has not been violated and the 
total number of particles is always conserved. Nevertheless, the continuity equation does 
not guarantee that the electronic density in each state does not go below zero or beyond 
two. 



As we mentioned in the introduction, the violation o 
observed for a different system in nuclear physics 21 



fermionic inequality has also been 



22j. In fact, there are earlier works in 



the classical BBGKY theory in which they studied the effect of nonlinearity introduced by 
truncation of the hierarchy, and showed the existence of instability in these coupled equations 
depending on the initial conditions of the system 
the classical collision integral can diverge 



42 



40 



41] . Other studies also indicated that 



43] . Such catastrophic behaviors of these 



coupled equations pose a valid question that why even highly advanced approximations based 
on the Green's function expansion not only fail to follow fundamental physical principles, but 
also lead to divergence, even though the total energy and number of particles are conserved. 

To analyze this phenomenon we center our attention to the basic properties of the BBGKY 
hierarchy and density matrices to find out how they are affected by different approximations. 
As we already showed, the employed approximations break the compatibility between Eq.flSJ 



and the approximated version of Eq. lflO]) and t he p artial trace relation 



Eq.(j5])) between F 



21] and Gherega et al. 



22J claimed this to 



and 7 does not hold any more. Schmitt et al. 
be the main source of the problem. 

On the other hand, it is obvious that the positive-semidefiniteness of density matrices 



1(3 



has also been violated. This problem may arise for one of the following reasons: 

i) If in Eq. ffTO]) the approximation functional of is built in a way that T does not 
necessarily stay positive-semidefmite, even though the initial 7 and T are positive- 
semidefmite (see cases C and D below). Therefore, regardless of whether the partial 
trace relation between 7 and T holds or not, there is no guarantee for 7 to be positive- 
semidefmite. 

ii) If the approximation functional of T^ 3 ^ is built in a way that the propagated T does stay 
positive-semidefmite (provided the initial 7 and T are positive-semidefmite), but since 
the Eqs. OH]) and ffTUl) are not compatible and the relation between T and 7 is ill-defined, 
the positive-semidefiniteness will not necessarily pass to the 7 (see case A below). 

It is not easy to impose positive-semidefiniteness on and even if it has such property, since 
its trace relation with 7 and T is broken this does not lead to the positive-semidefiniteness of 
7 and r. However, we are exploring variational approaches to impose this constraint during 
the time evolution, namely preventing the system not to be positive-semidefiniteness. 

Next, we use different test approximations to analyze the role of compatibility and 
positive-semidefiniteness in these unphysical results. 

A. Retaining Positive-semidefiniteness but not Compatibility. 

For the simplest case, 3b-CIF, the compatibility between second and first equation of the 
hierarchy is obviously lost. At the same time, it is easy to show that the time evolution of 
r is unitary in the sense that the solution to the equation of motion is of the form 

r(x 2 ,x' 2 ,t) = ^\ j $ j {x 2 ,t)$* j {x' 2 ,t). (31) 

j 

where $j satisfies the Schrodinger equation 

idt^j(X 2 ,t) = H 12 ^j(X 2 ,t) (32) 

Inserting Eq. (13 II) into Eq. ([()]) leads to time- independent geminal occupation numbers and 
provided the initial V is positive- semidefinite, it preserves this feature at all times (FigfS^a)). 
Note that this argument holds independent of possible time dependence in the two-particle 
Hamiltonian, H 12 . 
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In addition, substituting Eq. (l3Tj) in the right-hand side of Eq.(j9]) yields a simple linear 
equation for the time propagation of 7(xi,x' 1 ,£) which does not diverge. Nevertheless, 
because the partial trace relation (Eq.fll)])) between 7 and T is no more valid, natural orbital 



occupation numbers may not necessarily lie between zero and one [23( and as a consequence 
the density can gain negative values (Figfj^a)). With this example we come to the conclusion 
that positive-semidefiniteness of the T is not a sufficient condition for keeping the 7 positive- 
semidefinite. 



B. Retaining Compatibility and Positive-semidefiniteness. 

It is very insightful to look for an approximation that retains both of these properties 
and see how well that can perform in the calculation. The simplest possible approach is to 
replace the three-body collision integral term, S in Eq. lflO]) . with another term under the 
condition that the two equations become compatible. This condition, of course, does not 
determine a unique term, but to study its effect we consider the following simple expression 

S(X 2 , X' 2 , t) = (N- 2) (U(X 2 ) - U{X' 2 )) T(X 2 , X' 2 , t). (33) 

Thus, the second equation of the hierarchy becomes 

(idt - #12 + H V2 ,) T(X 2 ,X 2 ,t) = (N- 2)(U(X 2 ) - U{X' 2 )) T(X 2 ,X 2 ,t). (34) 

When in above equation, we put x 2 = x' 2 and integrate over x 2 we obtain the first equation 
of the BBGKY hierarchy, Eq.©. 

This approximation makes Eq.f JTQ|) linear. Moreover, it retains the compatibility between 
two equations, meaning the two approaches for solving the BBGKY hierarchy, discussed in 
section III Bl are equivalent. We call this the compatible approximation. The effect of this 
approximation is equivalent to magnifying the interparticle interaction by a factor of (iV— 1) 
in the 3b-CIF approximation. Hence, we can modify interparticle interaction accordingly in 
Eq.( l3~Tj) and use the same argument to prove that this approximation also keeps T positive- 
semidefinite. Now, as a result of compatibility, 7 = 7 r and this 7 immediately inherits 
positive-semidefiniteness of T and in this fashion, compatibility and positive-semidefiniteness 
are both incorporated. The results for this approximation are depicted in the Fig.(|3]) for a 
limited time scale. As expected, we do not observe any violation in the electronic density 
even after long propagation. 

18 



Compatible 

Exact 

TDHF 




Time [a.u] 

FIG. 3. Number of electrons in the leftmost site of a 4-site Hubbard model with compatible 
approximation. The exact and TDHF results are also given for comparison. Here, m, h are set to 
unity and Hubbard parameters are f7 = 0.1 and t = 1. The four electrons filled the two leftmost 
sites initially. 

Despite the well-behaved result, it should be mentioned that whereas the natural orbital 
occupation numbers cannot acquire negative values, they may still exceed one (the density 
at one point may develop beyond two particles). For example, we observed such violation 
for natural orbital occupation numbers in a tiny time region only when the initial state was 
far from equilibrium, i.e. four electrons occupying the two leftmost sites, and U was very 
large (U = 10). 



C. Retaining Compatibility but not Positive-semidefiniteness. 

In order to find out the role of positive-semidefiniteness and its importance, in the 
next step, we build an approximation that fulfills the compatibility but not positive- 
semidefiniteness. To accomplish this goal, we add a term, Z, with a partial trace summing 
up to zero, to the previous compatible approximation, Eq. fl33|) . Substituting this term in 
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Eg. (ITU]) we have 



(id t -H 12 + H ll2 >)T(X 2j X' 2j t) = 

(N - 2)(U(X 2 ) - U{X' 2 )) F(X 2 , X' 2 , t) + Z(X 2 , X' 2 , t) 



(35) 



where 



dx 2 Z(xi,x 2 ,x / 1 ,x 2 ,t) =0 



(36) 



which in most cases, Z is a nonlinear functional of density matrices. In this way, both 
levels of the hierarchy stay compatible but the additional term does not necessarily keep 
T positive-semidefinite. This term is not unique and in fact it can be derived for different 
approximations (see Appendix O for detailed derivation). 

Nonetheless, although these approximations retain the compatibility, they do not neces- 
sarily keep positive-semidefiniteness of RDMs and they eventually lead to divergence if we 
propagate them long enough in time. These results are the evidence that compatibility is 
not a sufficient condition to keep the equations bounded. 

D. Using Only the Second Equation (Eq. dlOj) ). 

As we discussed in section III Bl in details, there are two approaches to solve the BBGKY 
hierarchy in the second level. So far in all discussed approximations, we have propagated 
both Eqs.Q and ffTOl together; now, we turn to the second approach and solve only Eq. ffTOj) 
for above approximations. However, as we showed earlier, the total energy will not be 
necessarily conserved in this method unless the approximation makes the two equations 
compatible. Figure H] shows that the energy fluctuations in time are relatively large for all 
the approximations, which puts a question mark over the quality of this approach in general. 

Nonetheless in this way, the compatibility between equations is not an issue anymore 
which helps us to study the behavior of only the second equation. Our calculations show 
that, even though the fermionic inequality violation of 3b-CIF approximation has been cured 
in this way, the divergence in 3b-NIA or WC approximations has surfaced again and in fact, 
even the divergence time has not been improved in any of them. With these observations, 
we can confidently claim that the incompatibility is not the only cause of diverging behavior 
of the equations, and the role of positive-semidefiniteness should also be taken into account. 



20 



0.4 

d 



3b-CIF 

3b-NIA — - 
WC~ 



Exact 








20 



40 60 

Time [a.u] 



80 



100 



FIG. 4. Fluctuations of total energy of the 4-site Hubbard model for different approximations using 
only Eq, (|10j) . The exact result is also given for comparison. Here, m and h are set to unity and 
Hubbard parameters are U = 0.1 and t = 1. Four electrons filled the two leftmost sites initially. 

In these four subsections, we provided many test approximations fulfilling different con- 
strains to show that neither compatibility between equations nor positive-semidefiniteness 
of the approximations by itself, can keep the propagation of the RDMs inside the fermionic 
boundaries. In fact, although the nonlinearity introduced by most of the approximations to 
Eq. ljTU!) . might be the cause of divergence, the violation of fermionic inequality might exist 
even in the case of linear approximations as we saw in the 3b-CIF approximation. Therefore, 
it indeed takes both of these approximation constrains to tame such coupled equations. 

IV. SUMMARY AND CONCLUSIONS 

In this work, we pointed out the main challenges we face in order to decouple the hierarchy 
of equations for the time evolution of density matrices. In particular, we studied several 
approximations for the three-body RDM in terms of the one- and two-body RDMs. First, 
we showed that once an approximation for the three-body RDM is made, the equation of 
motion for the two-body RDM does no longer imply the validity of the first equation of 
the hierarchy. Therefore, in order to obey energy conservation it is necessary to solve the 
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equations of motion for the one-body and two-body density matrices simultaneously. 

Next, we studied numerical solutions for the fermionic Hubbard model for several of the 
decoupling schemes and compared them to exact results obtained from solving the time- 
dependent Schrodinger equation. We found that in many decoupling schemes the local 
electron density attains unphysical negative values in time and natural orbital occupation 
numbers in general do not remain between zero and one as is required for fermionic systems. 
Furthermore, in most of existing approximations the local electron density diverges, although 
total particle number and energy are perfectly conserved. We investigated whether this 
feature is cured by forcing the two lowest equations of the hierarchy to be compatible and 
found this not to be the case. 

We conclude that, a possible way to make progress in the application of the BBGKY hier- 
archy is to make sure that positive-semidefiniteness of the density matrices and the fermionic 
constraints on the occupation numbers is built into the equations. For example in the case 
of the time-dependent Hartree-Fock, there exists always a wave function corresponding to 
the first equation of the hierarchy and therefore the natural orbital occupation numbers can 
never be unphysical. This suggests a further study of the BBGKY equations in relation to 
(approximate) wave functions. Further progress can be made in linear response theory since 
in the linear response regime the nonlinearities are, by definition, removed and the BBGKY 
approach can be investigated further for the study of optical spectra. Work along these lines 
is in progress. 
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Appendix A: Energy Conservation Derivation 

Now we discuss in full details the energy conservation in the BBGKY hierarchy. According 
to Eqs. ffTTT) and (112"]) . the change in the total energy of the system, E(t) = Ei(t) + E 2 (t), is 
given by 

dE dE x dE 2 ,. . 
= t j £ (Al) 

dt dt dt K J 

where 

-j±= / dx' d t v(x!, t) 7 (x, x', t) \ x=x > + / dx' h(x!, t)dtj(x, x', t) | x=x , (A2) 

and 

dE 2 1 



dt dx l dx 2 U(X 2 )d t T(X 2 ,X 2l t). (A3) 

To proceed, we replace the time dervative of T from Eq. dlOj) . however for 7 in Eq. (lA2p . we 
can use either Eq.fjH]) or ( TT8|) . In our derivation, we use Eq. (fl8|) and the case of Eq.(j9]) is 
obtained as the special case by taking T = T at the end of the derivation. Let us start by 
evaluating (1A3[) . The equation of motion (fTOj) yields 

^ = ~J dx x diC2U(X 2 )(Vi + Vl- vl)r(x 2 ,x' 2: t)\ li2=v , 2/ 

= -— / dxidx 2 o r (xix 2 )(Vi - Vi,)r(xix 2 ,x / 1 x 2 ,t)|i=i' 
2% J 

= ^.J dxiViC/(x lX2 ) ■ [(Vi - ViOr(x 1 x 2 ,x' 1 x 2 ,t)| 1=1 ,] (A4) 

where we used the symmetry of V and performed a partial integration. Now, we consider the 
change in the one-body energy, Ei(t), primarily by evaluating the second term in Eq. (IA2j) . 
It reads 

J dx h(xf, t)d tlv (x, x', t) | x=x / = j dx (v(x, t)d t n T (x, t) + l -d t V ■ V 7r (x, x', t) | X=X ^A5) 
where in the second term of the right-hand side, we used 

j dxdx'5(x - x')V' 2 7 r (x, x', t) = -J dxdx'5(x - x')V • V' 7r (x, x', t). (A6) 
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On the other hand, from the equation of motion ([TBI) it follows 

id t Vi ■ Vr7r(xi,xi/,t)| 1=v = (Vi«(xi,*) ■ Vy - Vi>v(xi,t) ■ Vi)7 r (xi,x r ,t)|i =1 * (A7) 

1, 

2 



~Vi- [(V-V')(V 1 -Vi07 r (x,x / ,0| 1 , =1 ,] 



+(Vi-Vi/) y rfx 2 (f/(x 1 x 2 )-f/(x , 1 x 2 ))f(x 1 x 2 ,x / 1 x 2 ,t)| 1=r . 

The second term in the right-hand side of this expression is a total derivative and vanishes 
upon integration. Putting Eq.( 1A7j) back into ( 1A5jl and the resulted expression in Eq.( 1A2l) . 

we get 

dE f 

= / dx.(dtv(x,t)n r (x,t) +v(-x,t)(d t n r (-x,t) + V • j r (x,*))) 

y cZx 1 dx 2 (Vit/(x 1 x 2 ) ■ Vy - V 1 ^(x , 1 x 2 ) ■ Vi)r(x 1 x 2 ,x' 1 x 2 ,t)| 1=r . (A8) 

The second term under the first integral vanishes since the continuity equation holds. There- 
fore, for the total energy we finally arrive at 

dE f 

— = / dxd t v(x,t)n r (x,t) (A9) 
+ y dxidx 2 Vif/(xix 2 ) • [(Vi - Vi')(r(xix 2 ,x' 1 x 2 ,t) - f (xix 2 ,x' 1 x 2 ,t))|i=i'] . 

In the case that 7 evolves through Eq.Q, we must replace T with T and hence the final 
result reads 

dE f 

— = dx.dtv(x,t)n(x,t) (A10) 
which is exactly the energy conservation law as it is discussed in the main text. 

Appendix B: Computational Details for the Hubbard Model 

The real space description of Eqs. (jH]) and Eq. (fTUj) makes them very impractical for com- 
putational purposes. Therefore, they must be transformed to the matrix from, using an 
appropriate basis set. In the case of Hubbard model, the basis set is made of site-orbitals for 
each spin; for example, for a four-site Hubbard system, the basis set contains eight orbitals. 

Here we only concern with spin-compensated systems(z.e. S z — S = 0) , in which we can 



use many symmetries to eliminate the spins and simplify the equations of motion 44J. We 
define 
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7 (r 1 r' 1) t) = J2 7(ri^i, r^, t) (Bl) 

CTl 

and 

r(rir 2 , t) = ^ r(rie7ir 2 0- 2 , r'^r'^, t). (B2) 



where a denotes the spin coordinate. Therefore, the matrix from of Eqs.flHJ) and (flOl) for a 
M-site Hubbard model in site-orbital basis set read 

id^ij = t (7i (,-+!) (1 - 5ja/) + 7i(j-i)( 1 - fyi) - 7(i+i)i(l - <W) - 7(i-i)i(! - $1)) 

+ U (Tijjj -Fuji) (B3) 

and 

id t Tij k i = t(Tij^ + rji(l — S^m) + ^ij(k+i)i(^ — Ski) + ^ijk(i+i)(^ — <W) + ^ijk(i-i)(^ ~ $11) 

— r(i+i)jfcz(i ~~ <W) ~ r (i-i)jfcz(i — ^ii) — ri(j+i)fc«(i ~~ <W) - r^-i^^i — 5ji)) 
+ £f (r^fefc ^ — Tukk Sij) 

. ^ / (S)^^" 7 "^^" 7 " (3) CT i CT j CT « CT i CT j CT ™ (3) CT i <T 3 CT « CT i CT J CT " (3) CT i <T 3 CT ™ CT i CT J CT ™x , s 

+ ^ 2-~i v ijkklk + - 1 - ijJfcU ~ ijifcH ~~ ^ ijjklj J 

<Xj<Tj<Tn 

where indices denote the site numbers. The should be replaced by an approximation 
and then the spins can be integrated out using the same symmetry argument. These are 
the actual equations that we solve in this work. 

As it is mentioned in the main text, we choose a four-site Hubbard model where four 
electrons initially filled the two leftmost sites so that 

I 7m = 2, z = l,2 
7« = \ (B5) 
I Ji,j — 0, for the rest. 

Now, since this initial state is a Slater determinant formed by two site-orbitals, T has the 
exact form of Eq. ( 1231) and can be written as 

Fijki = li,k1j,i - ipi,ilj,k- (B6) 



For the initial state being the ground state of HF, we use Eq. (IB6p with the same argument 
and construct the V from 7. Evidently, we cannot use the same argument for any initial 
state. In the case of starting from the exact ground state, we extract the exact 7 and V and 
feed them into the equations. 
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Appendix C: Deriving The Zero- Traced Term 



Here, we derive the zero-traced expression, Z, which we introduced in section [III CI This 
additional term can be written for different approximations as 

Z{X 2 , X' 2 , t) = J dx 3 (t/(x lX3 ) + f/(x 2 x 3 ) - U(x[x 3 ) - f/(x' 2 x 3 )) r^( XlX2 x 3 , x / 1 x 2 x 3 , t) 
-F apP (^T) (CI) 

where the first term is the 3-body collision integral with I^ 3 ) substituted by an approximation 
and F app (7, T) is defined for that approximation as 

Tr X2 (F app ( 7 , T)) = J dx 3 rfx 2 (t/( Xl x 3 ) - f/(x' lX3 )) r^x^x,, x' lX2 x 3 , t). (C2) 

Evidently the partial trace of Z will be zero. 

As an example, we give the explicit form of F^-NiAd, T) as 

F 3b - NIA (l,r) = ([/(x lX2 ) - [/(x' lX2 )) {A( 7 (x 1 ,x / 1 ) 7 (x 2 ,x 2 ) - 7 (x 1 ,x 2 ) 7 (x 2 ,x / 1 )) 

+ x 2 )7(x 2 , x 'i) - #( x 2, x 2 )7( Xl , X i) 

+ 5 , (x 2 ,x' 1 )7( X i, X / 2 ) - 5f( X i, X ' 1 )7( X2 , X / 2 )} (C3) 

where 

5f(x / 1 ,x 1 ) = y rfx 2 7(x' 1 ,x 2 )7(x 2 ,x 1 ). (C4) 

In the same way, -F(7, T) can be obtained for other approximations but we do not present 
them here. Therefore, the full expression to replace the 3-body collision integral reads 

S(X 2 , X 2 , t) = (N- 2)(U(X 2 ) - U{X 2 )) r(X 2 , X' 2 , t) (C5) 
+ J rfx 3 (U (x lX3 ) + U (x 2 x 3 ) - U ( X ' lX3 ) - U ( X2 x 3 )) rgJ,( XlX2X3 , x' 1 x , 2 x 3 , i) 

-F app (7,r) 

where the first term is coming from the Eq. fl33|) . 
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